% y代表原序列
[y,Fs]=audioread('share.wav');
y=y(:,1); % 如果是双通道，可以使用该行代码，只考虑单通道
% 此处为了方便，将坐标的横轴显示成单位为s
T=1/Fs;
t=(1:length(y))*T;
subplot(3,2,1)
plot(t,y)
ylim([-2 2])
xlabel("时间/s");
ylabel("幅度");
title("原始语音信号时域图");
Y=fft(y);
subplot(3,2,2)
plot((1:length(Y))*2/length(Y),abs(Y))
xlabel("频率");
ylabel("幅度");
title("原始语音信号频域图");
ylim([0 1000])
% 高斯白噪声
noise=1*rand(size(y));
% 设计高通滤波器
Fs = 44100;  % Sampling Frequency

Fstop = 9600;            % Stopband Frequency
Fpass = 12000;           % Passband Frequency
Dstop = 0.0001;          % Stopband Attenuation
Dpass = 0.057501127785;  % Passband Ripple
dens  = 20;              % Density Factor

% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fstop, Fpass]/(Fs/2), [0 1], [Dstop, Dpass]);

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);

hpf=filter(Hd,noise);
noise_fft=fft(hpf);
% subplot(3,2,3)
% plot(t,hpf)
% subplot(3,2,4)
% plot((1:length(Y))*2/length(Y),abs(noise_fft));


y1=y+hpf;
%sound(y1,Fs)
% audiowrite("加噪.wav",y1,Fs);
subplot(3,2,3)
plot(t,y1)
% ylim([-0.6 0.6])
xlabel("时间/s");
ylabel("幅度");
title("加噪音语音信号时域图");
Y1=fft(y1);
subplot(3,2,4)
plot((1:length(Y1))*2/length(Y1),abs(Y1))
ylim([0 1000])
xlabel("频率");
ylabel("幅度");
title("加噪音语音信号频域图");
% 降噪处理
% 巴特沃斯低通滤波器
%通带截止频率
Wp=0.41;
% 阻带截止频率
Ws=0.5;
Rp=3;Rs=50;


[N,Wc]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(N,Wc,"low");


y2 = filter(b,a,y1);
% sound(y2, Fs);

% 去噪音语音信号时域图
subplot(3,2,5);
plot(t,y2);
ylim([-2 2]);
xlabel("时间/s");
ylabel("幅度");
title("巴特沃斯去噪音语音信号时域图");

Y2 = fft(y2);
% 去噪音语音信号频域图
subplot(3,2,6);
plot((1:length(Y2))*2/length(Y2),abs(Y2));
ylim([0 1000]);
xlabel("频率");
ylabel("幅度");
title("巴特沃斯去噪音语音信号频域图");

% 计算滤波前的信噪比
% 原始信号功率

psignal=P(y);
% 噪声功率
pnoise=P(y1-y);
% 信噪比
SNR_Original=10*log10(psignal/pnoise);
disp(SNR_Original)

% 计算滤波后的信噪比
% 降噪后的功率
psignal_after=P(y);
% 噪声功率
pnoise=P(y2-y);
SNR_After=10*log10(psignal_after/pnoise);
disp(SNR_After)



